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Abstract 

Background: Gliomas are the most common types of brain cancer, well known for 
their aggressive proliferation and the invasive behavior leading to a high mortality rate. 
Several mathematical models have been developed for identifying the interactions 
between glioma cells and tissue microenvironment, which play an important role in 
the mechanism of the tumor formation and progression. 

Methods: Building and expanding on existing approaches, this paper develops a 
continuous three-dimensional model of avascular glioma spatio-temporal evolution. 
The proposed spherical model incorporates the interactions between the populations 
of four different glioma cell phenotypes (proliferative, hypoxic, hypoglychemic and 
necrotic) and their tissue microenvironment, in order to investigate how they affect 
tumor growth and invasion in an isotropic and homogeneous medium. The model 
includes two key variables involved in the proliferation and invasion processes of cancer 
cells; i.e. the extracellular matrix and the matrix-degradative enzymes concentrations 
inside the tumor and its surroundings. Additionally, the proposed model focuses on 
innovative features, such as the separate and independent impact of two vital nutrients, 
namely oxygen and glucose, in tumor growth, leading to the formation of cell 
populations with different metabolic profiles. The model implementation takes under 
consideration the variations of particular factors, such as the local cell proliferation rate, 
the variable conversion rates of cells from one category to another and the nutrient- 
dependent thresholds of conversion. All model variables (cell densities, ingredients 
concentrations) are continuous and described by reaction-diffusion equations. 

Results: Several simulations were performed using combinations of growth and 
invasion rates, for different evolution times. The model results were evaluated by 
medical experts and validated on experimental glioma models available in the 
literature, revealing high agreement between simulated and experimental results. 

Conclusions: Based on the experimental validation, as well as the evaluation by 
clinical experts, the proposed model may provide an essential tool for the 
patient-specific simulation of different tumor evolution scenarios and reliable 
prognosis of glioma spatio-temporal progression. 

Keywords: Glioma, modelling, Spatio-temporal evolution, Cancer prediction, Tumor 
cells, Nutrients models, Proliferative, Hypoxic, Hypoglycemic, Necrotic regions 
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Introduction and background 

Cancer is the second most fatal disease worldwide after heart disease [1]. A can- 
cer cell evolves from normal due to genetic mutations, which abnormally alter the 
cell proliferation rate. In particular, glioma is a rapidly evolving type of brain can- 
cer, well known for its aggressive and diffusive behavior [2]. This diffusive invasion 
has lead several research efforts to explore the gliomas progression with the aid of 
mathematical diffusion equations [3-5], aiming to predict its spatial and temporal 
evolution. The high diffusion rate of glioma cells from the core tumor into the sur- 
rounding brain tissue often leads to treatment failure and tumor recurrence, even after 
the surgical resection. Gliomas vary from low- to high-grade, namely glioblastomas, 
which constitute the most malignant form of brain cancer, having an extremely poor 
prognosis. 

Several diverse factors influence glioma progression. Tumor cell migration in brain can 
be stimulated by the surrounding extracellular matrix (ECM) that is a complex mixture 
of macromolecules, some of which like the collagens play a structural role and others, 
such as laminin, fibronectin and vitronectin are important for cell adhesion, spreading 
and motility. ECM is remodeled by glioma cells that are able to degrade the surrounding 
tissue, through the secretion of matrix-degradative enzymes (MDEs), such as the plas- 
minogen activator and the family of matrix metalloproteinases. Essentially, all tumors go 
through an avascular stage of growth, where the nutrient (oxygen and glucose) supply is 
adequate. The nutrients are provided both by the intratumoral ECM and the surround- 
ings of tumor via diffusion. In a future growth stage, the intratumoral ECM is completely 
degraded, while the surrounding extracellular material is not close enough to supply the 
entire tumor body with the necessary nutrients via diffusion. Thus the tumor area is 
divided into layers consisting of different cell populations, namely proliferating, hypoxic, 
hypoglycemic and necrotic regions of tumor cells. Fully developed tumors contain a 
necrotic core with absence of nutrients, surrounded by hypoxic and hypoglycemic rings 
with limited access to nutrients and an outward proliferating cell ring having full nutri- 
ent supply. Since the tumor dependents on nutrients diffusion, its avascular growth is 
limited. Further growth and proliferation leads to its vascular stage, where cancer cells 
enhance the existing vascular network through angiogenesis. Furthermore, the metastatic 
final tumor stage succeeds the vascular stage. Since cancer development comprises these 
distinct stages (avascular, vascular and metastatic), the majority of research efforts are 
concentrated on exploring each one of these stages separately. In particular, identifying 
the avascular tumor dynamics constitutes the first crucial step towards fully vascularized 
tumors investigation [6] . 

In parallel to identification of tumor characteristics, the prediction of tumor growth 
and diffusion can lead to useful insight into the disease dynamics, which may improve 
clinical outcomes. To this respect, several mathematical and computational models have 
appeared in the literature, which investigate the mechanisms that govern gliomas pro- 
gression and invasion, with the aim of predicting its future spatial and temporal evolution, 
with or without the effects of therapy [7] . The models may constitute valuable tools for 
assisting the clinical practice towards the optimal individualized treatment, while facili- 
tating medical research analysis. Current modelling techniques focus on either simulating 
individual cell behavior, or modelling tumor of clinically significant size, exploring volume 
expansion [7,8]. Such models follow a discrete cell-based description that simulates cell 
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growth, or a continuum framework that deals with the evolution of local cellular densi- 
ties and is based on diffusion-reaction equations. Since both approaches reveal particular 
tumor aspects but also pose certain limitations, hybrid models have been developed 
to combine continuum and discrete variable descriptions, but they are still limited by 
model complexity and uncertainties in the mixing of models. Alternatively, the notion 
of multi-compartment continuum model appears to be effective in describing how sub- 
populations of various cell types proliferate and diffuse, while it remains computationally 
efficient [9]. 

Discrete models are able to incorporate cellular parameters and biological rules, as 
well as simulate the effects of chemotaxis, haptotaxis, cell to cell adhesion and address 
intracellular processes and intercellular communication issues. Cancer cells are usually 
represented by vectors of variables that determine their position and speed in space, 
their cell cycle phase, etc [10-13]. Moreover, these models are able to easily describe het- 
erogeneous cell populations. However, they are computationally intensive and fail to be 
efficiently initialized with tumors of clinically significant size that already consist of over 
1 billion cells. The discrete model of 1000 billion individual units is beyond the capability 
of most computers, even if simple interactions algorithms are exploited [14]. Thus, such 
models are limited to the exploration of small tumors, below the level of clinical detection. 

Continuum models represent alternative approaches in exploring tumor progression 
dynamics. Rather than considering each cell as a discrete individual, they assume tumors 
as spatial distributions of cell densities. As a result, computational requirements are sub- 
stantially reduced. If the population of tumor cells is obtained as a continuum variable, 
the system can be described at macroscopic level as a biological tissue sample. Continuum 
models are also capable of accounting for the dynamic changes of chemical ingredients 
(such as glucose, oxygen, MDEs, chemotherapeutic drugs, etc) in the brain tissue. Addi- 
tionally, they can take into account brain tissue heterogeneity and anisotropy. However, 
they lack a capability to simulate specific cellular factors, such as the genetic fingerprint of 
tumor cells and discontinuous changes, e.g. epithelial changes important for the invasion 
of tumor cells in adjacent tissue. Additionally, they are less efficient than discrete mod- 
els in describing cellular-scale flux factors such as chemotaxis, haptotaxis and cell to cell 
adhesion. 

Continuum approaches are most usually implemented by means of reaction-diffusion 
equations [15], where the cell proliferation factor is alternatively expressed by the 
Gompertz law, the second order polynomial equation, or the simple linear function 
[5,16]. Recent developments in mathematical modelling of gliomas are extensively ana- 
lyzed in [4]. First, a model of untreated glioma is described followed by the presentation 
of models including chemotherapy, or surgical resection. The modelling assumptions 
vary from relatively simple considerations dealing with homogeneous brain tissue, to 
those involving heterogeneous brain tissue with different diffusion rates of glioma cells 
in grey and white matter on a geometrically complex brain domain. In the approach 
of [16], a three-dimensional diffusion model has been developed to address brain het- 
erogeneity and anisotropy issues. The study introduces a mathematical framework 
incorporating alternative proliferation rate schemes and presents a solution of the dif- 
fusion equation through different numerical approximation of finite differences. Brain 
heterogeneity and anisotropy issues are dealt in the glioma growth model of [17], 
which is applied on real data. The proportions of white and grey matter, as well as 
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the diffusion tensors are extracted by normal brain atlases, avoiding DTI process- 
ing. The research effort in [18] supports the hypothesis that the change in pH of 
tumor microenvironment may provide an important mechanism for cancer invasion. 
Another reaction-diffusion model presented in [19], deals with the diffusive behavior 
of glioma growth. This homogeneous model includes the Neumann boundary con- 
ditions imposed by the skull on gliomas and especially on glioblastoma multiforme 
(GBM). In [20] growth and diffusion issues of U87 glioblastoma tumor spheroids 
are investigated. The model includes diffusion, proliferation and cell to cell adhe- 
sion parameters, which differ for each cell line. U87 glioma cells are also explored 
in [21], where a model that includes the role of chemotaxis, haptotaxis and adhe- 
sion is presented. According to this model cell migration is strongly depend on 
adhesion, haptotactic and chemotactic parameters, apart from the random diffusion 
factor. 

Expanding the potential of continuum mathematical models, some research efforts also 
engage the mechanical influence of tumor cells on the invaded tissue [5,22-26]. Many 
continuum approaches exploit the exponential Gompertz law [27] in order to model the 
proliferation of cancer cells. According to this law, the growth rates of cell populations 
are high at early stages of growth, while they are slower at later stages. The Gompertz 
tumor growth models take into account the fact that as the tumor size increases, the 
growth rate is reduced because of the reduction of available space and the limited supply 
of nutrients. The total cell population follows the Gompertz growth model in [28], where 
a two-compartment model of cancer cells population dynamics is proposed. The model 
consists of a set of ordinary differential equations, which simulate proliferating and qui- 
escent cell populations including the transition rates between them. Moreover, Gompertz 
law and Gompertzian parameters available in the literature are exploited in the models of 
[29] and [30]. 

Recently, hybrid models have been developed to overcome particular limitations of 
distinct approaches. These models exploit the continuum method in order to simulate 
the tumor microenvironment, while cell to cell interactions are implemented using the 
discrete approach [14,31-37]. Some of these models are multiscale, incorporating the 
molecular and cellular level [6,38-40]. However, hybrid models are still constrained by 
the discrete-model limitations in which they do not consider tissue-level associations and 
they still have a significant computational cost. 

Multi-compartment continuum models form alternative continuum approaches that 
are expected to overcome some of the conventional modelling limitations and are 
able to include multiphase heterogeneous populations and heterogeneous extracel- 
lular matrix [41]. According to this model philosophy, cells are grouped based 
on phenotype depending on their access to the necessary nutrients, or the level 
of differentiation [42]. Such approaches are computationally efficient and effective 
in describing how subpopulations of various types of cells proliferate and diffuse. 
Invasion, proliferation, changes in phenotype and necrosis are readily expressed as 
additions or subtractions from the cell densities within each compartment. In [3] 
a continuum multi-compartment model is proposed investigating how tumor cells 
interact with their tissue microenvironment. The dispersal and interactive popu- 
lation changes of normoxic and hypoxic glioma cells, vascular endothelial cells, 
diffusible angiogenic factors and necrotic cells are implemented using diffusion 
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equations. Furthermore, the study demonstrates how different proliferation and 
diffusion rates of glioma cells result to increasing degrees of mitoses, hypoxia- 
induced neoangiogenesis and necrosis, features that characterize the different glioma 
grades. 

In this paper, a new continuum three-dimensional spherical and multi-compartmental 
mathematical model of avascular glioma growth in an isotropic and homogeneous 
medium is proposed. This model expands the multi-compartmental approach of [3] by 
integrating a new cell compartment that represents the hypoglycemic cell population. 
The initial model is further extended through the exploitation of the microenvironmental 
variables of the hybrid approach in [32], where apart from oxygen, we introduce a sec- 
ond vital nutrient, (i.e. glucose). The proposed model consists of heterogeneous tumor 
cell populations incorporating the interactions between four different glioma cell phe- 
notypes into distinct cellular compartments, namely proliferative, hypoxic, hypoglycemic 
and necrotic, as well as their tissue microenvironment. Moreover, the model includes 
the effect of the host tissue, i.e. the extracellular matrix (ECM), the matrix-degradative 
enzymes (MDEs) and nutrients (oxygen and glucose) concentration on tumor cell prolifer- 
ation and invasion, through different microenvironmental compartments. The proposed 
approach employs continuum variables (cell densities, concentrations of chemical ingre- 
dients) governed by diffusion principles. It incorporates the effects of concentration 
for both important nutrients in the tumor microenvironment, namely oxygen and glu- 
cose. Notice that the majority of existing approaches deal with a single nutrient, which 
most often stands for oxygen [14,21,24,26,32,35,37,40,43] and causes the development 
of hypoxic cell regions. Furthermore, existing approaches employing both vital nutri- 
ents (e.g. [1,6,34,39,44]) consider only the combined effect of their availability, which 
accounts for the generation of only quiescent (non-proliferating) and necrotic cells. The 
novelty of our proposed model lies in the fact that apart from the combined effect, it 
takes into consideration the separate and independent impact of each nutrient in tumor 
growth, which leads to the formation of a new cell population (apart from hypoxic, 
or quiescent and necrotic), namely hypoglycemic cells, with a different metabolic pro- 
file. Moreover, this model includes innovative variation of particular factors, such as the 
local cell proliferation rate, the variable conversion rates of cells from one category to 
another and the nutrient-dependent thresholds of conversion. The model is validated 
through a comparison with specific experimental results of glioma models available 
in the literature and additional evaluation by clinical experts. The results demonstrate 
the models efficiency, providing an essential tool for the patient-specific simulation of 
different tumor evolution scenarios and reliable prognosis of glioma spatio-temporal 
progression. 

The paper is organized as follows. In Section "Mathematical modelling", the differ- 
ent cell and chemical compartments are described in detail as building blocks of the 
model, coupled with the appropriate differential equations. Section "Results and model 
validation" presents and discusses experimental results of tumor growth for different 
combinations of diffusion and proliferation rates and for various evolution times. More- 
over, model validation is performed through comparisons with experimental results of 
glioma models available in the literature, as demonstrated in the same section. Finally, sig- 
nificant conclusions are drawn in Section "Discussion and conclusions", where potential 
future work aiming to model improvement and completion is also discussed. 
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Mathematical modelling 

Current glioma modelling techniques, either apply discrete cell-based tumor growth sim- 
ulations focusing in individual cell behavior, or follow continuum approaches dealing with 
the evolution of tumor cell densities. The discrete models are computationally feasible 
only if the initial number of tumor cells is relatively small, which is not possible for a tumor 
of clinically significant size that is detectable in current imaging techniques (MRIs, CTs, 
etc). Taking this into consideration, the rationale behind choosing to develop a continuum 
glioma evolution model in this paper relates to the need of dealing with virtual tumors 
of clinically visible size and density, consisting of large numbers of cancer cells. This 
model is initialized by a virtual spherical glioma tumor and simulates multiple cellular- 
microenvironmental interactions aiming to predict its three-dimensional spatio-temporal 
evolution. 

The evolution of glioma can be described by the following reaction diffusion equation 
[4,5]. 

9 C 

— = V-J(C,t)+S(C,t)-T(C,t) (1) 
dt 

where / represents the flux/motility of tumor cells, S stands for their net proliferation and 
T is the factor concerning the death of cancer cells due to microenvironmental conditions 
(e.g. absence of nutrients) and the beneficial effect of the applied treatment. 

The motility of the tumor cells can be expressed as the result of four different fluxes 
/ = J Random + JHaptotaxis + Jchemotaxis + J J Adhesion* namely the random cell diffusion, as 
well as the diffusion due to haptotaxis, chemotaxis and cell to cell adhesion [20,21]. 
J Random = DrVC concerns the random cell diffusion, which is constant in a homoge- 
neous medium and JHaptotaxis = Dj-iapt ' CVf stands for the haptotaxis that is the directed 
migratory response of cells to gradients of non-diffusible chemicals, such as the extra- 
cellular matrix. Moreover, Jchemotaxis = ^chem ' CVGl is the tendency of cancer cells to 
move in the direction of the nutrients (glucose) gradient and J Adhesion = C • K(C) cor- 
responds to the movement due to cell to cell adhesion, according to which cells adhere 
to each other when they are close, while push apart when they are too compressed. In 
these equations Dr, Dnapu ^chem> are the random diffusion, the haptotactic and the 
chemotactic coefficient respectively, C, f, Gl the corresponding cancer cells, extracel- 
lular matrix and glucose concentrations and K represents the adhesive force of tumor 
cells. 

The net proliferation rate of tumor cells S can be described by either a simple linear 
function S(C, t) = p • C, the second order polynomial equation S(C, t) = p • C • (1 — -£-), 
or the Gompertz law«S(C, t) = p-C'ln(^) [5,16], where C is the concentration of cancer 
cells, p corresponds to the proliferation rate constant and C m stands for the tissue cell 
maximum carrying capacity. 

Cancer cells are capable of degrading the surrounding tissue, i.e. the extracellular matrix 
(ECM), which constitutes a crucial part of the invasive process. This degradation is 
achieved by means of the matrix-degradation enzymes (MDEs) that are produced by the 
cells and governed by the diffusion-reaction equation [32]. 

dm 9 , x 

— = D m V 2 m + ii-C-X-m (2) 
dt 

where m represents the concentration of MDEs, D m is the MDEs diffusion coefficient, C 
stands for the concentration of cancer cells, fi is the MDEs production rate and X stands 
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for the MDEs natural decay rate. The concentration change of ECM is simulated by the 
following equation [32] . 



where e represents the concentration of ECM and 8 is the ECM degradation rate. 

The tumor cells in order to survive, proliferate and diffuse need nutrients including 
oxygen and glucose. Many of them follow the anaerobic metabolism (glycolysis) instead 
of the normal aerobic metabolism that utilizes oxygen, even under of adequate oxygen 
presence [6,45]. The cells that follow the glycolytic pathway do not use oxygen; however 
they need a much larger quantity of glucose, in order to produce energy, in comparison to 
those that have an aerobic metabolism [34]. Hence, the tumor requires both oxygen and 
glucose, which are provided by the blood vessels and diffuse to the tumor area through the 
intratumoral and surrounding ECM. The nutrient production, diffusion and consumption 
are expressed by the diffusion-reaction equation [32] . 



where n represents the concentration of the vital nutrient, D n is the nutrient diffusion 
coefficient, /3 is the nutrient production rate through the ECM, y stands for the nutrient 
consumption rate by the cancer cells and a is its natural decay rate. 

Based on these equations, the proposed model expands the idea of glioma cell diffu- 
sion and proliferation taking into consideration different cell categories according to their 
metabolic profile, as well as the critical effect of two vital nutrients, which are consumed 
by the tumor cells. 

Proposed model description: concepts and assumptions 

The proposed continuum mathematical approach models the three-dimensional spher- 
ical spatio-temporal evolution of an avascular glioma in an isotropic and homogeneous 
medium. The model incorporates specific interactions between the heterogeneous pop- 
ulations of cancer cells and their microenvironment, which affect tumor growth and 
invasion. It takes into consideration the separate and independent effect of the concen- 
tration changes for both of the important vital nutrients in the tumor microenvironment, 
namely oxygen and glucose, leading to the formation of distinct cell populations with 
different metabolic profile and invasion properties. This innovative feature of the pro- 
posed cancer model can deal with different cell-profile changes for each nutrient, thus, 
allowing the formation of different areas in the cancer neighborhood. In other related 
approaches, when the concentration of one of the two nutrients is sufficiently reduced 
(no matter which) only one new cell population (e.g. quiescent cells) appears before 
leading to necrosis. The proposed model consists of continuous variables, described by 
diffusion equations. The invasion of tumor cells, which is governed by the diffusion 
principles is based on the following equation: 



where C is the concentration of cancer cells at tumor point (r, 0,0) (spherical coordi- 
nates) at time £, D represents the coefficient of random cell diffusion and/ stands for 
the cell net proliferation rate. The proliferation rate in this model is based on the second 
order polynomial equation, since it takes in account the ratio of the local cell density and 




(3) 



an 9 

— = D n V 2 n + 0 • e- y • C -a • n 



(4) 



9C(r, 0,0,0 
dt 



= DV 2 C(r, 0, 0, t) +/(C(r, 0, 0, t)) 



(5) 
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the maximum carrying capacity of the tissue, in order to limit the cell proliferation. Thus, 
when the tissue reaches its maximum carrying capacity, the tumor cells in this region, 
stop proliferating and remain quiescent, which is reasonable due to the lack of space. Fur- 
thermore, our model allows for variable local-cell proliferation rate due to the nutrients 
availability. 



In order to decrease the complexity of the model, we express the invasion of tumor cells 
only due to random diffusion, i.e. haptotaxis, chemotaxis and cell to cell adhesion have 
been neglected. The latter mechanism, according to which cells exert adhesive forces to 
each other, is observed only on the cellular scale. Since our model deals with the macro- 
scopic consideration of cell densities, it is reasonable to ignore these forces. Moreover 
chemotaxis, namely the tendency of cancer cells to move in the direction of the nutri- 
ents, has been neglected in the current model. This consideration does not significantly 
degrade the approximation of real tumor cell-motility, since the tendency of the cells to 
move in the direction of nutrients supply is implicitly taken into account in our model. 
More specifically, due to random diffusion, the tumor cells tend to move from the areas of 
large cell concentrations towards less crowded regions, which however hold higher nutri- 
ent concentrations because of their lower consumption rate in comparison to the regions 
of large cell populations; hence, tumor cells are indirectly forced move towards nutri- 
ents. Similarly, the effects of haptotaxis, which is the directed migratory response of cells 
towards non-diffusible chemicals, e.g. ECM, are implicitly considered; through random 
diffusion, cells tend to move towards regions of lower densities, which happen to hold 
higher concentrations of ECM. 

The proposed continuum model of spherical shape is multi-compartmental. The evolu- 
tion of tumor mass is modeled by means of coupled diffusion equations, which simulate 
the different cell densities. The compartmental approach is exploited to track glioma cell 
subpopulations based on viability and phenotype. Thus, this model includes a prolifera- 
tive region, along with a hypoxic, hypoglycemic and a necrotic one. The cellular size of 
all phenotypes is assumed to be the same and the cell diffusion coefficients are consid- 
ered constant, since the tumor is assumed to evolve inside an isotropic and homogeneous 
medium. Additionally, coupled continuum components are developed in order to account 
for the changes in the local environment that affect cells within each compartment. These 
components concern the diffusion and consumption of vital nutrients, concentration and 
degradation of the host tissue, i.e. ECM, as well as production and diffusion of MDEs. 
The diffusion coefficients of nutrients and MDEs are assumed constant within all tumor 
heterogeneous regions and outside the tumor as well. The incorporated nutrients that 
determine tumor cell viability and affect proliferation are the oxygen and glucose. 

Oxygen and glucose are transported to the tumor region from blood vessels and diffuse 
through the ECM located within the tumor and its surroundings, which is destroyed by 
MDEs produced by the tumor cells as the tumor grows. Initially, the size and cell den- 
sity of the tumor are small, thus all cells are sufficiently supplied with oxygen and glucose 
through the intratumoral and surrounding extracellular material [1,6,21,26,31,36]. At this 
time all cancer cells are normoxic and proliferative, since the diffusion of nutrients is ade- 
quate for their survival and proliferation and most of them follow the aerobic metabolism. 



f(C(r,0,4>,t)) = p.C(r,e,4>,t).Q- 




(6) 
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Furthermore, these cells compete for space with other neighboring cells; if the local tis- 
sue reaches its maximum carrying capacity, they stop proliferating and remain quiescent, 
while additionally they reduce their metabolism [32]. 

As time proceeds, the tumor radius increases and the concentrations of oxygen and glu- 
cose decrease in the central part of the tumor, closer to its core. At a certain time, each 
nutrient reaches a critical concentration value, below which it is not sufficient to meet 
the needs of all cells. Nutrients do not simultaneously fall below their critical concentra- 
tions, but are depleted independently depending on their consumption rates. Once the 
first nutrient, say oxygen, attains its critical concentration, the local tumor cells become 
hypoxic at a specific rate; a hypoxic region appears in the tumors center. Similarly, when 
the concentration of the other nutrient, i.e. glucose, falls below a critical level insufficient 
to meet energy needs of all cells, local cells turn to hypoglycemic and a hypoglycemic 
region appears in the central part of the tumor. The hypoxic and hypoglycemic cell pop- 
ulations neither proliferate nor die, while their metabolism changes. Specifically, the 
hypoxic cells perform increased glycolysis [45], which requires much less oxygen; the 
amount of oxygen consumed in hypoxic regions is not precisely known and could actually 
reach five times less the oxygen consumption in proliferative cells [32,46] . However, since 
anaerobic metabolism is far less efficient in producing energy compared to the normal 
aerobic one, hypoxic cells eventually consume over 10 times more glucose for maintaining 
the same energy turn-over [34]. As far as hypoglycemic cells are concerned, they reduce 
their metabolism, i.e. consumption of nutrients [6,32,35] and remain quiescent (they stop 
proliferating). Moreover, hypoxic cells are more migratory than normoxic proliferative 
cells [42], thus their diffusion coefficient is higher while they secrete larger amounts of 
degrading enzymes. When hypoxic cells migrate to an area with adequate oxygen, they 
convert back to normoxic at a specific constant rate; exactly the same applies to the case 
of the hypoglycemic cells, with respect to glucose [3]. The cells that remain in a state of 
oxygen and glucose adequacy in a hypoxic or a hypoglycemic area continue proliferating; 
however, as time proceeds these cell populations decrease, since they increasingly convert 
to hypoxic or hypoglycemic respectively. Eventually, in the hypoxic and/or hypoglycemic 
areas the total proliferation rate is reduced to the appropriate levels. As the tumor evolves, 
both hypoxic and hypoglycemic populations grow in the central part of the tumor at the 
expense of proliferative cells, to finally form a hypoxic and/or a hypoglycemic core within 
the tumor. 

When the vast majority of cells have turned to hypoxic or hypoglycemic in the cen- 
tral part of the tumor, one of the nutrients (oxygen or glucose) tends to vanish and 
tumor cells start dying, i.e. become necrotic. Thus a necrotic region appears close to the 
tumor core, where nutrients can no longer diffuse due to the long distance from the tumor 
surroundings. Additionally, necrotic cells appear in a tumor area where both nutrients 
have achieved their critical concentrations, regardless of hypoxic and/or hypoglycemic 
cell proportion. Actually, any of the nutrients can be the first to extinct at the tumor cen- 
ter depending on the initial concentrations, their diffusion coefficients and consumption 
rates. The necrotic cells do not have a metabolism, do not proliferate and do not move, 
since they are not alive. Moreover, necrosis can also result from the viable proliferative 
cells due to their contact with the necrotic ones, at a particular rate [3]. In addition to 
necrosis, the natural cell death process is also present in the form of, apoptosis [1]. How- 
ever, apoptosis has been neglected in this model as to keep the complexity to a manageable 
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level, since the apoptosis rate is low enough in comparison to the necrosis rate and does 
not significantly affect the cell proliferation and invasion processes. 

Model equations for tumor-cell densities 

Based on the above considerations of tumor progression, our model eventually considers 
tumor volume formed in four areas as in Figure 1, composed of the necrotic core, the 
hypoxic along with hypoglycemic rings around it, and the outer proliferating ring. Each 
of these areas consists mostly of cells belonging to the corresponding category. As long as 
the tumor radius increases, the hypoxic, hypoglycemic and necrotic radii increase as well. 
The tumor reaches a steady state when the hypoxic, hypoglycemic and necrotic regions 
occupy most of the entire tumor body. 

In the model of Figure 1, four distinct glioma cellular compartments are developed. 
Each compartment constitutes a continuum equation, coupled with the remaining cellular 
compartments via common factors between them. The invasion of proliferative cancer 
cells follows the diffusion-reaction equation. The cell proliferation rate is expressed by the 
second order polynomial function, which takes into consideration the tissue maximum 
carrying capacity and diminishes the rate inversely according to the local cell density raise. 
As mentioned before, in a region of maximum cell density, the tumor cells apart from the 
fact that they do not proliferate, they reduce their metabolism (i.e. nutrient consumption) 
as well. Although these quiescent cells have temporarily stopped proliferating, they are 
included in the population of proliferative cells, since they still stay on a normoxic and 
normoglycemic state just as proliferative cells. The cell density equations are based on 




Figure 1 Tumor regions in avascular evolution. Tumor regions in avascular evolution, where P is the 
Proliferative, H is the Hypoxic, Q is the Hypoglycemic and N is the Necrotic zone. 
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[3] and have been expanded to include the new cell compartment that represents the 
hypoglycemic population, as well as the effect of oxygen and glucose concentration in cell 
type (proliferative, hypoxic, hypoglycemic, necrotic) and proliferation rate. The equations 
governing the growth and invasion of proliferative cells are presented below. 

J = V(D C • (1 - T)VC(r,0,4>,t)) +/(C(r, 0,0,0) 
at 

+ H(r, 0, 0, t) • g h • (1 - n h ) + Q(r, 0, 0, t) • g q • (1 - ^) (y) 

- C(r, 0, 0, t) -b h - n h - C(r, 0, 0, ^ 

- C(r, 0, 0, £) • Af(r, 0, 0, 0 • * w 

This equation involves several terms describing the various cellular processes 
discussed above, where C, H, Q, N denote the concentrations of proliferative 
(normoxic/normoglycemic), hypoxic, hypoglycemic and necrotic cancer cells respectively 
at tumor point (r, 6, 0), at time t, D c represents the coefficient of random proliferative cell 
diffusion and T stands for the fraction of the local cell density to the maximum carrying 
capacity, given by: 

= C(r, 0, 0, t) + H(r, 0, 0, t) + Q(r, 0, 0, t) + N(r, 0, 0, t) 

— p \ ) 

Furthermore, bh, b q are the conversion rates of proliferative cells to hypoxic and hypo- 
glycemic respectively, which unlike [3] are variable proportional to the oxygen or glucose 
relative concentrations. Additionally, a n stands for the constant of conversion rate of pro- 
liferative cells (as well as hypoxic and hypoglycemic) to necrotic due to the contact with 
the necrotic region and g^, g q , are the constant conversion rates of hypoxic and hypo- 
glycemic cells to proliferative, respectively. Moreover, n^gl q represent the thresholds for 
the proliferative cells to turn to hypoxic or hypoglycemic and in this model are related 
to the critical nutrients concentrations. The function/ in Equation 7 stands for the local 
proliferation rate of tumor cells and corresponds to an innovative net proliferation rate 
variation, since in our model depends on the local oxygen and glucose availability. It can 
be modeled as: 

/(C(r, 0, 0, t)) = O con • Gl con • p • C(r, 0, 0, t) • (1 - T) (9) 

where p corresponds to the proliferation rate constant and O con , Gl con are coefficients that 
reduce the proliferation rate in relation to the local oxygen (n) and glucose (gl) concentra- 
tion respectively and are given as in the following equations. The concentrations n andg/ 
have been calculated including the current consumption by the cells at each time point, 
while ntmp and gkmp are the corresponding concentrations before nutrient consumption. 
Additionally, no and glo represent the maximum concentrations of the nutrients which 
initially exist in the tumor surroundings. In particular, 

^ i (u , fo, w(r,0,<M) > 0 n tmp(r,0,<l>,t) 

O con = l-(b h +a n )-n h , n h = \ and b h = (1 )/20 

[1, n(r,0,<l>,t) < 0 n 0 

(10) 



Glcon =l-(b q +a n )-gl q , gl q = 



0, g i {r ,e,<p,t) > o W) )/20 

1, gl(r,6,<p,t) <0 q glo 

(11) 
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The invasion of hypoxic cells is governed by diffusion, following the equation below: 

dH(r,0,4>,t) 



= V(D h • (1 - r)V/f(r,6>,0,O) 

(12) 



dt 

+ C(r,6>,0,£) • b h • n h - H(r,0,<l> 9 t) • g h • (1 - n h ) 
- H(r, 0, 0, i) -a h ' n hn - H(r, 0, 0, t) • a gih • gl hn 
-H(r,0,<l>,t)>N(r,0,<l>,t)'a n 



where represents the coefficient of random hypoxic cell diffusion, is the variable 
(proportional to b n ) conversion rate of hypoxic cells to necrotic due to the lack of oxygen, 
while cigih, is the conversion rate of hypoxic cells to necrotic, due to the absence of glu- 
cose, rih n , glh n represent the thresholds for the hypoxic cells to turn to necrotic due to the 
absence of oxygen, or glucose respectively, modeled as: 



a H(r,0,4>,t) < 0Q 

U ' C(r,0,cf>,t)+H(r,0,(f>,t)+Q(r,0,(f>,t) U ' y j j _ 

-, H{r,6,(t),t) nq ana S l hn — 

A ' C(r,0,</>,t)+H(r,6,<f>,t)+Q(r,0,(/>,t) > U ^ 



0, n h = 1, gl q = 0 

h n h = 1, gl q = 1 



Similarly the diffusion of hypoglycemic cells is expressed by the following equation: 

8Q(r,0,d>,t) 

d J =V(D q • (1 - DVQ(r, 0,0,0) 

+ C(r, 0, 0, £) • V gl q - Q(r, 0, 0, t) • g q . (1 - gl q ) (13) 
- Q(r, 0, 0, t)>a q > n qn - Q(r, 0, 0, 0 • a giq • g/^ 
-Q(r,0,<l>,t)>N(r,0,<l>,t).a n 

where D q represents the coefficient of random hypoglycemic cell diffusion, a q is the con- 
version rate of hypoglycemic cells to necrotic due to the lack of oxygen, while a g \ q is the 
variable (proportional to b q ) conversion rate of hypoglycemic cells to necrotic, due to the 
absence of glucose. Furthermore, n qni gl qn represent the thresholds for the hypoglycemic 
cells to turn to necrotic due to the absence of oxygen or glucose respectively, described in 
the form: 



glqn = 



°> C(r^,0,O+2(r^l0,f)+Q(r,^,0,O - 075 » _ { °> g l q ~ "h - 0 

1 QW4,t) Q75 ana n in - . _ 

A ' C(r,0,</>,t)+H(r,0,4>,t)+Q(r,0,<f>,t) > "' /D [ L > &q ~ L > n h — 1 



Finally, necrotic cells do not diffuse and thus their concentration changes follow the 
rather simple equation: 

dN(r,0,<l>,t) 

= C(r, 0, 0, t) • N(r, 0, 0, t) • a n 

dt 

+ H (r, 0, 0, f) • flfc • + H (r, 0, 0, f) • ^ • ^ ^ 
+ /f (r, 0, 0, 0 • N(r, 0, 0, 0 • a n + Q(r, 0, 0, f) • ^ • n qn 
+ Q(r, 0, 0, t) • ^ • ^ + Q(r, 0, 0, £) • N(r, 0, 0, f) • a n 

The factors participating in Equation 14 have all been mentioned above. 
Model equations for nutrient concentrations 

The concentration-changes of oxygen and glucose are similarly governed by diffusion 
equations, which includes their production through the ECM, their natural decay and 
their consumption by the tumor cells. The oxygen concentration equation, as well as the 
new glucose concentration compartment is based on [32]. Both of the nutrients equations 
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have been expanded to involve the different consumption rates by each of the distinct cell 
types (proliferative, hypoxic and hypoglycemic). More specifically: 

dn(r,0,d>,t) 9 

— ^— = D n V 2 n(r, 0, 0, t) + p„ • e(r, 0, 0, f)-a n > n(r, 0, 0, t) 

ot 

- Yen • C(r, 0, 0, t) - y hn • H(r, 0, 0, t) - y qn • Q(r, 0, 0, t) (15) 
3rf(r, 0, 0, t) 9 

- d J = D gl V 2 gl(r, 0, 0, t) + ^ • e(r, 0, 0, 0 

- • gKr, 0, 0, 0 - y cg i • C(r, 6>, 0, t) 

- • H(r, 0, 0, 0 - y qg i • Q(r, 0, 0, t) (16) 

where ^/ represent the concentration of the oxygen, glucose at tumor point (r, 0, 0), at 
time £, Dg/ are the respective diffusion coefficients, fi n , fi g i are the nutrients produc- 
tion rates through ECM, a n , a g i are the nutrients natural decay rates, y cn , y cg i stand for the 
nutrients consumption rates by the proliferative cancer cells, y/^, yhgl are tne nutrients 
consumption rates by the hypoxic cells and respectively y qn , y qg i are the nutrients con- 
sumption rates by the hypoglycemic cells. In the tumor area where local cell-density has 
reached the tissue maximum carrying capacity, the normoxic/normoglycemic cells have 
been considered in this model to consume half of the nutrients' amounts [32]. 

The concentration change of ECM because of its degradation by MDEs produced by 
the tumor cells, is demonstrated by the following equation (based on [32]): 

d jMlM = _ 8 • m(r, 0, 0, t) • e(r, 0, 0, t) (17) 
dt 

where e is concentration of ECM at tumor point (r, 0,0), at time t, m represents the 
concentration of MDEs and 8 is the ECM degradation rate. MDEs are subjected to the 
principle of diffusion, produced by proliferative, hypoxic and hypoglycemic cells and are 
characterized by a natural decay rate. The MDEs equation is based on [32] and has also 
been expanded to include the distinct tumor cell categories: 
dm(r, 0,0,0 9 

v = D m V 2 m(r, 0, 0, t) + /x c . C(r, 0, 0, t) 

dt (18) 

+ fi h • H(r, 0, 0, t) + fi q - Q(r, 0, 0, t) - X • m(r, 0, 0, f) 

where D m is the diffusion coefficient of MDEs, fi Cf /jl^ fi q are the production rates by 
the proliferative hypoxic and hypoglycemic cells respectively and X represents the MDEs 
natural decay rate. 

Parameter estimation 

One of the major limitations of mathematical models concerns the parameter initializa- 
tion using literature-based reference values. These values are often difficult to obtain, 
since most of them are patient specific and vary significantly among clinical and exper- 
imental cases. However, the existing wide range of values for each parameter allows 
models to incorporate different abnormal cases and examine many alternative outcomes. 

A continuum glioma model can be initialized based either on the actual tumor geome- 
try derived from medical images, or on a virtual spherical tumor [1], which approximates 
a real one, in terms of its imaging detectable size. In imaging modalities such as CT 
and MRI, the tumor is detectable above a certain size and density [3,4]. According to 
the literature, the density threshold associated with the detectable boundary of a glioma 
in a medical image has been estimated to C2D = 400ce//s / mm 2 , which leads to C^d = 
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(V^2d) 3 = 8 * 10 3 cells I mm 3 , while the smallest detectable size corresponds to a radius 
of r = 1 cm. According to medical experts tumors of radius smaller than 1 cm (e.g., 
r = 5 mm) can be diagnosed. Yet, taking into consideration that oxygenated cells are 
limited to a distance of less than 0.1 — 0.2 mm from blood vessels, most of the initially 
detectable tumors (r > 5 mm) should already contain at least a hypoxic, even a necrotic 
region. Thus, in order to address tumorogenesis at early stages, our model is initialized 
with a virtual tumor of a smaller radius (i.e. r = 1 mm), under the assumption that the 
initial tumor consists only of one cell category, namely proliferating cells without a signif- 
icant obliquity. Accordingly, the initial tumor density has been considered to be slightly 
higher than the corresponding minimum detectable threshold, i.e. Co = 10 4 cells /mm 3 . 
This constitutes the initial condition for the current model, which involves two bound- 
aries; the inner one (r = rn) corresponding to the tumor border and the outer boundary 
(r = r\) taken far enough from the tumor and defining the area of the medium to be 
invaded by the tumor. The initial condition for tumor cell density is described as: 



The invasion rate of glioma varies among patients depending on the grade of the dis- 
ease. A wide range of invasion rates has been observed in clinical experience from serial 
MRIs and can be found in the literature in the form of several diffusion coefficient val- 
ues. Similarly, for another patient-specific parameter namely the constant of tumor cell 
proliferation rate, a wide range of values has been reported. Different combinations of the 
diffusion and proliferation rates determine the different glioma grades. For instance, con- 
cerning a four grade malignant glioma (glioblastoma), a primary one is characterized by 
low diffusion and high proliferation, while a glioblastoma progressing from lower grade is 
determined by higher diffusion but lower proliferation rate [3]. In this paper, various com- 
binations of diffusion and proliferation rates have been considered, in order to explore 
different invasion and growth cases associated to various glioma grades. As far as the dif- 
fusion of hypoxic cells is concerned, since they are more migratory than normoxic cells 
[42], their invasion rate is higher; here their diffusion coefficient is taken one order of 
magnitude higher compared to the normoxic cells value. On the other hand, due to the 
lack of related data, the diffusion coefficient of hypoglycemic cells has been considered 
the same as the one of proliferative cells. 

Initially, apart from outside the tumor, ECM is located within the tumor area as if it 
has not been completely degraded by MDEs. In our consideration, the initial intratumoral 
ECM concentration is considered one order of magnitude lower than the ECM maxi- 
mum concentration at the tumor surroundings. As the tumor grows, intratumoral ECM is 
destroyed by MDEs, which are produced by the proliferative, hypoxic and hypoglycemic 
cells at specific rates. Since hypoxic cells are more invasive than normoxic cells, their 
MDEs secretion is considered larger and this is reflected in our model by accounting a 
twice higher MDEs production rate. Regarding hypoglycemic cells, they are associated to 
the same MDEs production rate as proliferative cells. The initial concentration of MDEs 
has been considered zero outside the tumor boundary and nonzero inside the tumor. 
More specifically, since the initial tumor density is lower than the tissue maximum car- 
rying capacity, the initial intratumoral MDEs concentration is considered lower than the 




(19) 
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maximum reference value, equal to the ratio of the tumor density over the maximum 
carrying capacity. 

The initial tumor is sufficiently provided with nutrients through the remaining intratu- 
moral ECM, as well as from the surrounding ECM via diffusion. The maximum nutrients 
concentration exists in the area outside the tumor, which serves as the nutrient source. 
In the current model, nutrient concentrations are reduced towards the central part of 
the tumor [1] according to a normal distribution, depending on the distance from tumor 
periphery as demonstrated by the following equation: 

Kotumon gktumor = Oo> gk) ' e~ dr (20) 

where no, glo represent the initial oxygen and glucose concentrations outside the tumor 
and d r is the distance of the current point from tumor periphery. The diffusion, as well 
as production, consumption and natural decay-rate parameters of the different chem- 
icals (i.e. nutrients, ECM and MDEs) are taken from related research studies, such as 
[1,3,21,32,33,35,46-49]. Particularly, regarding the metabolic rates of the various cell 
phenotypes, hypoglycemic cells are assumed to consume nutrients at half the rate of pro- 
liferating cells. On the contrary, since the vast majority of hypoxic cells are considered to 
follow the glycolytic pathway, they are assumed to consume 10 times more glucose, but 
1/5 less oxygen [42], compared to normoxic cells, the majority of which (approximately 
60%) are considered to follow the aerobic metabolism in this model. Some parameters 
were hard to obtain and they were estimated according to related known values. The 
parameter values that were used in the equations of our model are summarized in Table 
Al contained in an additional file (please see Additional file 1). 

Boundary conditions for all tumor cell densities are determined by the first set of the 
following equation pairs. The same boundary conditions apply for MDEs concentration. 
Finally, nutrients and ECM are determined by the second equation pair. 



dr ~ U ' 1 ~ U 

C(r,0,4>,t) = 0, r = n 



dn{r,6,(t),t) _ n _ n 

dr r -" (21) 

n(r,6,(j),t) = n 0 , r = n 



Results and model validation 

Simulation results 

The compartmental equations of the proposed model along with its solutions have 
been numerically approximated by means of a finite differences method, namely For- 
ward Euler. Several simulations have been performed using combinations of diffu- 
sion and proliferation rate constants, for the various evolution times of 1, 2, 3, 6, 
9, 12, 15 and 18 months. More specifically, the results concern the combinations 
of [low diffusion, low proliferation], [high diffusion, low proliferation], [low diffu- 
sion, high proliferation], [high diffusion, high proliferation] and [medium diffusion, 
medium proliferation]. The model is initialized with an early-stage tumor of radius 
ro = 1 mm, smaller than the minimum detectable size consisting only of prolifer- 
ating cells, without any other assumptions on the initial radius of the various cell- 
population regions. Furthermore, the quiescent cells due to the lack of space have been 
included in the population of proliferative cells, since they still remain on normoxic 
and normoglycemic conditions. The simulation results of cell and nutrient densities 
with respect to the distance from the center of the tumor, for the combination of 
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low diffusion (D = 0.005 mm 2 /day) and high proliferation rate (p = 0,04/ day), 
corresponding to a primary high grade glioma, i.e. a glioblastoma from its first clinical 
detectability, are graphically illustrated in Figures 2, 3, 4 and 5 for the evolution times 
of 1, 6, 12 and 18 months. More results concerning the intermediate evolution times of 
2, 3, 9 and 15 months are illustrated in an additional file (please see Additional file 2). 
This additional file also contains the graphical evolution of the low diffusion-high prolif- 
eration glioma evolution in the early stage of 1 day after its assumed detection, where all 
cancer cells remain proliferative since no hypoxia and/or hypoglycemia have appeared yet. 
Finally, all results concerning radius of regions for the various cell categories for the dif- 
ferent times stages and the diffusion-proliferation pairs of interest are depicted in Table 1. 

Figure 2 demonstrates the evolution of the primary glioblastoma 1 month after its 
assumed detection, when its radius has increased from r = 1 mm to r = 3.2 mm. It 
can be noticed that due to the high proliferation rate and the low invasion, cell den- 
sity rapidly increases in the tumor center, which leads to the appearance of hypoxia 
(blue curve and region in Figures 2(a), 2(b)). The hypoxia is also verified in Figure 2(c) 
reflecting the oxygen inadequacy in the central part of the tumor. Concerning the other 
diffusion-proliferation combinations demonstrated in Table 1, the cells remain prolifer- 
ative throughout the tumor, implying that there is still adequacy of oxygen and glucose. 
This proves that as D gets larger, the tumor is dominated by its diffuse extent that reduces 
cell density in its central part, thus it takes longer for the proliferative cells to fall into 
hypoxic and hypoglycemic conditions. 
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Figure 2 Model simulation of a low diffusion-high proliferation tumor for 30 days (1 month). Model 


simulation of a low diffusion-high proliferation tumor for 30 days (1 month): (a) cell density with respect to 


the tumor radius, (b) tumor cells dispersal in a 2D section, (c) oxygen and (d) glucose concentrations with 


respect to the tumor radius. 
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r(mm) r(mm) 

(c) (d) 

Figure 3 Model simulation of a low diffusion-high proliferation tumor for 1 80 days (6 months). Model 
simulation of a low diffusion-high proliferation tumor for 1 80 days (6 months): (a) cell density with respect to 
the tumor radius, (b) tumor cells dispersal in a 2D section, (c) oxygen and (d) glucose concentrations with 
respect to the tumor radius. 



Simulation results of 6 months after the assumed tumor diagnosis (low diffusion, high 
proliferation) in Figure 3, demonstrate that the hypoxic cell density at the tumor cen- 
ter has increased. Additionally, a hypoglycemic cell population has appeared (magenta 
curve in Figure 3(a)) due to the glucose inadequacy, which is shown in Figure 3(d). Glu- 
cose concentration is further decreased since hypoxic cells consume much more glucose 
in order to survive. Apart from hypoxic and hypoglycemic cell populations, a necrotic 
core of approximately 3 mm, has been formed (red curve and region in Figures 3(a), 
3(b)). This is followed by the hypoxic and hypoglycemic regions around it and the exter- 
nal proliferative region. Here the hypoglycemic zone is incorporated inside the hypoxic 
band. The hypoxic, hypoglycemic and necrotic populations have increased at the expense 
of proliferating cells (green curve in Figure 3(a)), which now constitute the minority for 
the tumor radius of r < 5 mm. Moreover, at the same time hypoxic, hypoglycemic and 
necrotic regions have appeared for the high diffusion, high proliferation case (Table 1), 
where the tumor radius already exceeds the value of r = 4 cm. This case could represent 
a progressive glioblastoma. 

After 12 months (Figure 4), the hypoxic and hypoglycemic areas have grown while the 
separate zones have become more distinct, since the cells that compose each zone pre- 
vail on all other cells. It can be clearly observed that proliferating cells no longer exist in 
the central part of the tumor. Additionally, it is shown that the hypoglycemic region is 
still incorporated within the hypoxic one. Moreover, the glioma represented by medium 
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Figure 4 Model simulation of a low diffusion-high proliferation tumor for 360 days (1 2 months). 

Model simulation of a low diffusion-high proliferation tumor for 360 days (1 2 months): (a) cell density with 
respect to the tumor radius, (b) tumor cells dispersal in a 2D section, (c) oxygen and (d) glucose 
concentrations with respect to the tumor radius. 



diffusion, medium proliferation rates has progressed to grade IV, containing a necrotic 
central region, representing a secondary glioblastoma. 

Finally, Figure 5 illustrates tumor growth after 18 months, where all the distinct regions 
have further increased in size and density from inside to outside. The formation of the 
extended necrotic core, surrounded by the hypoglycemic and hypoxic regions and the 
outward proliferating zone are clearly depicted. 

Table 1 and Figure 6 demonstrate the relation between the diffusion coefficient and the 
tumor-growth rate. More specifically, it can be observed that in all cases the expansion 
rate is temporarily reduced after the first 3 months, while later it is increased proportion- 
ally to the diffusion-coefficient value concerning the medium and high diffusion cases, to 
finally reach a steady growth state. Regarding the low diffusion cases, the tumor growth 
rate is stabilized after its first reduction. 

Model validation 

Model validation constitutes an essential part of mathematical modelling and it can be 
performed using either medical data stemming from clinical cases, or experimental model 
results available in the literature. The model in this paper is based on a virtual spherical 
tumor, so that its validation is more appropriate using experimental model results derived 
from similar existing research efforts. However, the same model can be initialized from a 
real tumor anatomy derived from MRI imaging. Nevertheless, the comparison with actual 



Papadogiorgaki etal. Theoretical Biology and Medical Modelling 2013, 10:47 
http://www.tbiomed.eom/content/1 0/1 /47 



Page 19 of 31 




clinical data requires the use of nonhomogeneous media considerations, which is outside 
the scope of this paper. The results of the current study are compared to the correspond- 
ing outputs of the continuum model in [3] and the hybrid model in [1]. In both cases, high 
correlation was revealed between simulated and experimental results. 

First, the experimental results drawn from graphs presented in paper [3] are compared 
to simulations results of the proposed model. Grade association with size and time is 
depicted in Figure 7 for two different diffusion-proliferation combinations in comparison 
to respective results of [3]. Specifically, Figure 7(a) demonstrates a primary glioblastoma 
(grade IV from its first detectability) and corresponds to low diffusion and high prolif- 
eration rate, while Figure 7(b) simulates a secondary glioblastoma, (progressing to grade 
IV from lower grade) and stands for medium diffusion and medium proliferation rate. 
For both cases, the initial tumor has a density equal to Co = 10 4 cells /mm 3 and a radius 
of r = 1 mm (as considered in all the above-presented simulation results). This plot 
demonstrates the high agreement between the results of the proposed model and the 
experimental results of [3] concerning spatio-temporal evolution and grade association 
of glioma, even though the proposed model simulates an avascular tumor. 

Table 2 depicts the time from the typical detectable size of r = 1 cm to the typical fatal 
size of r = 4 cm for all the five different combinations that were investigated in this paper. 
Additionally, the time intervals for the same D/p pairs are shown in this table for model in 
[3], which were approximated from the color-scale graph of predicted survival times for a 
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Table 1 Model simulation for 1, 2, 3, 6, 9, 1 2, 1 5 and 18 months 



Model simulation results 
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0.005 


0.3 
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0 
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0 


0 


20 
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20 



Radius (mm) of each cell population (Proliferative: P r , Hypoxic: H r , Hypoglycemic: Q r , Necrotic: N r ) with respect to each 
diffusion-proliferation combination. 



spectrum of in silico patients. It can be observed that almost in all cases, the time that the 
tumor needs to reach the fatal radius of 4 cm according to the proposed model predic- 
tion is within the time-intervals of model [3] for the corresponding diffusion-proliferation 
tumor types. 

In the model of [3] the simulated T2-detected tumor corresponds to a total cell density 
greater than the T2-threshold of detection, which according to this paper is approximated 
to 16% of the maximal tissue-carrying capacity. Following this rationale for comparison 
purposes, our model is also initialized with a spherical tumor of a density lower than 
the maximum carrying capacity, yet higher than the minimum detectable value of [3] 
(0.16 • C m < Co < C m ), i.e. Co = 10 5 cells / mnv > . Figures 8(a), 8(c) and 8(e) illustrate the 
simulations of each grade of glioma, represented as a plot of cell density for the different 
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Figure 6 Glioma grade as a function of tumor size and time. Tumor grade as a function of tumor size and 
time for each of the five diffusion-proliferation (D - p) combinations: (a) low diffusion-low proliferation, (b) 
high diffusion-low proliferation, (c) low diffusion-high proliferation, (d) high diffusion-high proliferation, (e) 
medium diffusion-medium proliferation. Blue part: grade II radius, green part: grade III radius, red part: Grade 
IV radius. 



cell populations (normoxic, hypoxic and necrotic), with respect to the distance from the 
center of the tumor. In Figure 8(e) the hypoglycemic cell population has been included 
in the normoxic cell population. Comparable curves derived from [3] corresponding to 
gliomas of similar grades and sizes are shown in Figures 8(b), 8(d) and 8(f). These figures 
demonstrate that a similar growth pattern is followed in both models, with respect to the 
different tumor regions in each glioma grade. 

The deviations regarding cell-densities between the two models signify their differ- 
ences and potential limitations. In particular, the proposed model, lacks a consideration of 
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Figure 7 Glioma grade as a function of tumor size and time in comparison to [3] results. Tumor grade 


as a function of tumor size and time in comparison to results of [3] for two different diffusion-proliferation 


pairs, (a) corresponds to a primary glioblastoma produced by D = 0.01 mm 2 /day, p = 0.033/day, (b) 




illustrates a secondary (progressive from lower grade) glioblastoma produced by D = 0.06 mm 2 /day, 




p = 0.01 3/ day. The dotted curves correspond to comparative results of [3]. Yellow part: minimum T2 




detectable radius, blue part: grade II radius, green part: grade III radius, red part: grade IV radius. 





angiogenesis, but considers the impact of glucose gradual reduction on the cell prolifera- 
tion rate. Moreover, the model of [3] is initialized with a tumor of a tenfold radius, namely 
the smallest T2 detectable size of r = 1 cm, compared to the proposed model. This large 
initial size with high cell-density within its entire volume significantly differentiates the 
final cell densities, i.e. when the tumor has reached its fatal size. 

In order to validate the proposed model for an initial tumor of a small radius r < 
1 mm, the results are compared to respective outcomes of model [1], which simulates 
avascular spherical tumor growth based on cell proliferation that depends on the con- 
centrations of oxygen and glucose. Model results are demonstrated for a growing tumor 
far below its clinical detectable size. The initial concentration values of both nutri- 
ents outside the tumor are the same used in our model (no = 0.28 • 10~ 9 mol/mm 3 , 
glo = 16.5 • 10~ 9 mol/mm 3 ). Since the initial tumor radius here is considered very small, 
r = 0.05 mm, the concentrations of both nutrient inside the tumor are assumed to have 
their maximum values, equal to the reference concentration. The model in [1] does not 
include a diffusion factor for the tumor expansion, but its radius increases corresponding 
to cell proliferation. For comparison purposes, the diffusion coefficient value of our pro- 
posed model has been considered to be proportional to the mean area that each tumor cell 
occupies. The cell proliferation rate, simulates the tumor expansion due to the increasing 
space of proliferating cells. Additionally, the proliferation-rate constant in this case has 
been taken from [1] as p = 0.12/ 'day. These modified diffusion and proliferation rates and 
the tumor initialization for r = 0.05 mm, lead to a similar tumor growth as the one in [1]. 



Table 2 Tumor evolution time interval from 1 to 4 cm 



Tumor evolution time 


D(mm 2 /day) 


0.005 


0.3 


0.005 


0.3 


0.04 


p(fday) 


0.0025 


0.0025 


0.04 


0.04 


0.015 


Modeling] 


6-9 years 


3-18 months 


2-3 years 


3-6 months 


1-2 years 


Proposed Model 


10 years 


1 7 months 


2 years 


5 months 


17 months 



Time interval for tumor evolution from the typical detectable size of 1 cm to the typical fatal size of 4 cm. 
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Figure 8 Glioma grades represented as a plot of cell density for the different cell populations. 

Simulations of each grade of glioma represented as a plot of cell density for the different cell population 
(normoxic, hypoxic and necrotic) with respect to the distance from the center of the tumor. Green curve: 
normoxic cells, blue curve: hypoxic cells, red curve: necrotic cells, (a), (c) and (e): proposed model produced 
by different diffusion-proliferation pairs, (b), (d) and (f): comparable results of [3]. 



The results are presented in Figure 9, where the tumor volume expansion for 28 days in 
comparison to the respective results of [1] are depicted. Concerning this growth, the radii 
of quiescent and necrotic regions in respect to the tumor radius are illustrated in Figure 10 
and compared to the experimental results of [1]. The quiescent radius corresponds to 
the hypoxic along with the hypoglycemic region of the proposed model Figure 10 clearly 
demonstrates a very good agreement between the experimental data and the simulated 
growth curves, resulting from our model prediction. 
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Figure 9 Tumor volume expansion in respect to time. Tumor volume expansion in respect to time. Green 
curve: experimental results of [1], blue curve: proposed model. 



Discussion and conclusions 

Gliomas form an important class of brain cancer with high mortality rate. Mathemati- 
cal models are increasingly used to understand and predict their behaviour. Identification 
of the avascular tumor dynamics constitutes the first crucial step towards the investi- 
gation of fully vascularized tumors. However, using current modelling techniques one 
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Figure 1 0 Quiescent and necrotic radius in respect to different tumour radius. Quiescent and necrotic 
radius in respect to different tumour radius. Blue curve: quiescent radius of proposed model, red curve: 
necrotic radius of proposed model, blue and red dashed curves: corresponding experimental results of [1]. 
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must choose between simulating individual cell behavior and modelling tumors of clin- 
ically significant size. Both the continuum and discrete tumor models pose particular 
limitations. In order to overcome significant restrictions, a multi compartment contin- 
uum model appears to be effective in describing how subpopulations of various types of 
cells proliferate and diffuse, while it is computationally efficient. To this aim, a new con- 
tinuum three-dimensional spherical mathematical model of avascular glioma growth in 
an isotropic and homogeneous medium is developed and proposed in this study. This 
model simulates malignant cell proliferation and invasion behavior, incorporating the 
interactions between the heterogeneous populations of four different glioma cell pheno- 
types, namely proliferative, hypoxic, hypoglycemic and necrotic, as well as their tissue 
microenvironment. Our model is capable of capturing intercellular interactions, spa- 
tial cell population heterogeneity and phenotype differentiation. It is validated through 
multiple simulations in order to determine sensitivity to changes in important model 
parameters, specifically, the proliferation and cell migration rates [4]. 

The different cell populations are implemented by means of distinct cellular com- 
partments based on [3], which has been expanded towards the integration of tumor 
microenvironment with the inclusion of chemical compartments. Unlike the majority of 
related existing approaches that involve a single nutrient (oxygen) to guide tumor growth, 
our model incorporates the effects of the concentration changes of two nutrients (oxygen 
and glucose). The novelty of the proposed model lies in the fact that apart from the com- 
bined effect, each nutrient separately affects tumor cell proliferation and viability, while a 
new cell population is formed due to the lack of glucose, characterizing hypoglycemic cells 
by a different metabolic profile than hypoxic ones. Moreover, the effects of ECM, MDEs 
and the two vital nutrients, on cell survival, proliferation and invasion are simulated based 
on the equations of [32] (for oxygen, ECM and MDEs), which have been extended to 
include the different production and consumption rates corresponding to the distinct cell 
types. In essence, our model builds on relevant existing models by combining their capa- 
bilities and expanding their exposure and applicability to realistic glioma scenarios. All 
model variables, e.g. cell densities and concentrations of chemical ingredients are of con- 
tinuous form described by diffusion-reaction equations. Several simulations have been 
performed using various diffusion and proliferation rate combinations corresponding to 
glioma reference values. The model results are presented for different evolution times in 
order to express the evolutionary characteristics of tumors. Furthermore, the model is val- 
idated through comparisons with specific experimental results of glioma models available 
in the literature, revealing high concordance of our model results with different aspects 
of glioma spatio-temporal evolution and tumor initialization phases. 

Since the clinical cases of substantial significance are above the tumor size of 0.5 — 1 cm 
radius, our model is initialized with a much smaller tumor-radius r = 1 mm contain- 
ing only normoxic proliferative cells, as to consider the evolution of tumor composition. 
Our model predicts tumor radius and composition after its clinical diagnosis, while it also 
allows tumor initialization at a size far below the detectable in MRI images, so that it can 
approximate the actual tumor cell-densities since its infancy. As it was shown, the model 
enables identifying the size of the different cell regions, such as proliferative, hypoxic, 
hypoglycemic and necrotic, as well as their effect on the overall tumor growth. Therefore, 
the proposed model may prove useful for determining particular tumor growth param- 
eters, which are hard to obtain from patient data but are crucial in the development of 



Papadogiorgaki etal. Theoretical Biology and Medical Modelling 2013, 10:47 
http://www.tbiomed.eom/content/1 0/1 /47 



Page 26 of 31 



therapeutic strategies for cancer treatment. Moreover, the validation results demonstrate 
that our model can work effectively and provide the tumor-growth prediction initialized 
either from a very small or a larger tumor size. This enhances the applicability of the 
model irrespective of the stage of diagnosis. 

In this study the effect of hypoxia as a hallmark of aggressive tumor behavior often 
met in glioblastomas that has been engaged in the model through nutrients availabil- 
ity/consumption has been associated with cell proliferation and invasion. According to 
medical experts this is of great significance in clinical practice, since it is associated with 
resistance to therapy, poorer survival and more malignant tumor phenotypes [2]. Thus, 
one important utility of models is to simulate a quantitative link between tumor growth 
kinetics and the hypoxic burden of the tumor. To this aim, our mathematical model 
describes the spatial and temporal evolution of glioma in terms of concentration of malig- 
nant tumor cells. A major strength of this formalism is its potential prospective nature. 
Although it does not include the prediction of the tumor recurrence location, our method 
is potentially applicable to orient patient-specific definition of glioma margins. 

Furthermore, since our model is able to detect cell densities of even below lcell/rnrn 3 , 
it can effectively characterize the border line regions, i.e. the semi-cancerous regions, 
as well as the rest extra-tumoral suspicious area, which potentially consist of healthy 
tissue along with scattered cancer cells. In addition, the area that is occupied by an 
extra-tumoral edema is likewise considered as suspicious area by the radiologists, since it 
contains high concentrations of matrix-degradative enzymes (MDEs), which are known 
to degrade the extracellular material (ECM) and provide a favourable environment to can- 
cer cells in order to migrate. Our model incorporates the variables of both the MDEs 
and ECM, thus it is able to detect those extra-tumoral suspicious areas containing high 
MDEs-concentrations and relatively low ECM-concentrations. 

Despite its utility and effectiveness, the proposed model in its present form bears cer- 
tain limitations. Tumor cells with various phenotypes initially develop asymmetric tumor 
morphologies, since they have different invasion rates and they migrate within an inho- 
mogeneous medium; thus they are non-uniformly distributed. However, they eventually 
form a circular shape so that the spherical shape of the initial tumor used in our model 
constitutes a successive real tumor approximation. Furthermore, this model is homoge- 
neous, i.e. it does not take into consideration brain heterogeneity regarding tumor cell 
invasion. As mentioned in [5], there are several indications that glioma diffuses faster in 
white than in grey matter. The cell diffusion coefficient considered in the current model 
is constant throughout the entire medium, while it should vary depending on the invaded 
tissue area. Along the same direction, model initialization is performed with a virtual 
tumor that lacks the specific spatial localization of a real tumor surrounded by particular 
brain anatomic structures (i.e. white and grey matter), where tumor cells would have dif- 
ferent invasion rates. Thus, one limitation of the current model is the lack of support for 
medium heterogeneity, which however can be considered at a further stage of develop- 
ment. Moreover, a real tumor consists of at least three distinct regions of different texture 
properties (e.g. necrosis, proliferative, edema), where chemicals have different diffusion. 
The model presented does not consider variable diffusion properties for the nutrients due 
to intratumoral heterogeneity, since diffusion coefficient values are very difficult to obtain 
and additionally, such a differentiation would not offer a significant accession to model 
results. 
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Cell to ECM interaction (i.e. haptotaxis) and cell to cell adhesion play their specific 
roles in tumor growth and morphological dynamics. Cell to cell adhesion is important at 
the early stages of glioma development, while at later stages tumor growth is dominated 
by cell to ECM and nutrients interactions. Cell to cell adhesion could be incorporated 
in a macroscopic continuum approach, but has not been included in our model, since 
it concerns the cellular scale that is not explored in this work. Regarding haptotaxis, it 
has also been neglected as a distinct diffusion term, but it has been implicitly taken into 
account through random diffusion of tumor cells towards areas of lower cell concen- 
tration, where ECM concentration is higher. Similarly, chemotaxis has been implicitly 
incorporated through the factor of random diffusion. As another remark, the proposed 
model does not evaluate the impact of gene expression changes and mutations in cellu- 
lar behavior, nor does it incorporate the effects of the epidermal growth factor receptor 
(EGFR) and the associated molecular pathways on tumor growth dynamics. 

Overall, the major effort in this work is to demonstrate the results of the changes in 
the dominant cellular kinetics, i.e. diffusion and proliferation, considering the interac- 
tions between tumor cells and their microenvironment. Hence, the current effort forms 
a comprehensive modelling of the change in the tumor cellular proliferative and invasive 
phenotype due to the dynamic cell interactions with the microenvironment that contains 
the vital nutrients. As such, it may miss several key factors that determine cell phenotype, 
which nevertheless can be readily explored in future models. 

Future efforts include the incorporation of several parameters aiming to overcome the 
limitations. A straight forward direction for model amendment is to address brain tis- 
sue heterogeneity and anisotropy issues, in order to take into account the different cell 
diffusion rates, depending on the invaded anatomic structure (i.e. white or grey matter). 
This can be easily performed by adopting the topography of imaging modalities (i.e. MRI) 
similar to other related studies [17]. In the same direction we should consider the model 
initialization by a real tumor geometry extracted from corresponding medical images. 
Additionally, primary tumor formation in terms of proliferative, hypoxic and necrotic 
regions, derived from imaging texture, or biopsy data could also be integrated in this 
initialization. This primary intratumoral heterogeneity will allow the integration of vari- 
able nutrient diffusion properties. Subsequently, the implementation of additional flux 
factors, such as diffusion due to haptotaxis, chemotaxis and cell to cell adhesion should 
be addressed. The current model considers the effect of nutrient-concentration changes 
but it does not take into consideration the increase of waste concentrations, such as 
lactate and H + ion, which act as cell proliferation inhibitors and threaten cell survival. 
The incorporation of those inhibitors is also considered essential in future developments. 
Furthermore, the effects of EGFR concentration and related pathways on tumor growth 
dynamics at a multicellular level will constitute a significant addition. Since the vascular 
morphology of gliomas has led to the hypothesis that the formation of new blood vessels 
is essential for tumor growth, future plans will target the investigation of angiogenesis 
within the proposed model. 

According to the medical experts in our team, this model is successful in capturing avas- 
cular tumor growth observed clinically. It allows performing patient-specific simulation of 
different tumor evolution scenarios, towards reliable prognosis of glioma spatio-temporal 
progression. Since it explicitly incorporates the microenvironment interactions of prolif- 
erative, hypoxic and hypoglycemic glioma cells, along with the formation of necrosis, it 
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allows evaluating clinically significant tumor sizes, while it is efficient in describing the 
dynamics of glioma tumors visualized with medical imaging. This mathematical model 
provides a means to simulate tumor development scenarios, which may lead to a bet- 
ter understanding of how altering fundamental parameters can influence brain tumor 
progression; thus it may constitute an important research tool for clinical assessment 
[50]. In this regard, it can be exploited in cooperation with in-vivo and in-vitro mod- 
els after establishing the initial diagnosis via a tumor biopsy or surgery. In-vitro invasion 
assays are important tools for investigating the tumor-matrix interactions and the effects 
of extracellular macromolecules on these interactions [51-53]. Since such studies are 
carried out through tightly-controlled experimental conditions, they allow the explicit 
determination of primary variables, e.g. nutrients initial concentrations [54,55]. On the 
other hand, in vivo implanted xenographs (human brain tissue) in animal model varia- 
tions can effectively investigate the reproducible cell migration and the tumor invasion 
into living non-neoplasmatic brain regions [56]. Thus, both in vitro and in vivo stud- 
ies of human tumor parts allow the examination of tumor behavior deriving important 
patient specific parameters, such as the invasion, the proliferation and the nutrients 
consumption rates [57]. The proposed model can then be initialized, based on these 
parameters in order to attain a detailed and accurate short- and long-term patient specific 
prediction. 

An important area of model exploration concerns the changes of tumor microenvi- 
ronment in the context of a specific treatment. The treatment of tumors in the central 
nervous system represents a formidable challenge, further magnified by the fact that 
the brain is isolated by the blood brain barrier, rendering the delivery of high doses of 
chemotherapeutic agents, or gene vectors a very difficult task [58]. This leads to unac- 
ceptable toxic systemic levels of drugs. The approaches proposing combinations of local 
tumor delivery, within the tumor mass and/or surrounding cavity, such as immune- 
stimulation and cytotoxic gene therapies, appear very promising in achieving the effective 
treatment of glioblastoma [59]. Additionally, the introduction of new technologies, such 
as the microchip and convection enhanced drug delivery, will enable the local delivery 
of treatments, such as drugs and therapeutic gene vectors within the tumor mass, as 
well as the surrounding area, where the infiltrating tumor cells are localized [60,61]. This 
provides the motivation for simulating and modelling the targeted cancer treatment sce- 
narios. The multi-compartmental continuous nature of our model enables the integration 
of therapy-related issues in the form of either a simple factor, or a pharmakokinetic- 
pharmakodynamic model, which will investigate the beneficial impact of treatment in 
disease progression. 
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